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Abstract 



Markov chain Monte Carlo (MCMC) simulations are commonly employed for es- 
timating features of a target distribution, particularly for Bayesian inference. A fun- 
damental challenge is determining when these simulations should stop. We consider a 
sequential stopping rule that terminates the simulation when the width of a confidence 
interval is sufficiently small relative to the size of the target parameter. Specifically, 
we propose relative magnitude and relative standard deviation stopping rules in the 
context of MCMC. In each setting, we develop sufficient conditions for asymptotic va- 
lidity, that is conditions to ensure the simulation will terminate with probability one 
and the resulting confidence intervals will have the proper coverage probability. Our 
results are applicable in a wide variety of MCMC estimation settings, such as expecta- 
tion, quantile, or simultaneous multivariate estimation. Finally, we investigate the finite 
sample properties through a variety of examples and provide some recommendations to 
practitioners. 

Keywords. Batch means, Bayesian computation, fixed-width confidence intervals, 
sequential estimation, sequential stopping rules, strong consistency. 

1 Introduction 

Markov chain Monte Carlo (MCMC) methods allow exploration of intractable probability 
distributions by constructing a Markov chain whose stationary distribution equals the de- 
sired distribution. A major challenge for practitioners is determining how long to run an 
MCMC simulation. Many experiments employ a fixed-time rule to terminate the simula- 
tion; that is, the procedure terminates after n iterations, where n is determined heuristically. 
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Indeed, some simulations are so complex that this is the only practical approach, but that 
is not so for most experiments. 

Alternatively, many practitioners use convergence diagnostics to determine if n is suffi- 



ciently large (for a review see Cowles and Carlin, 1996). Although practical, these methods 
are mute about the quality of the resulting estimates (Flegal et al. 2008). Moreover, they 



can introduce bias directly in to the estimates (Cowles et al. 1999). 

We instead advocate terminating the simulation when an estimate is sufficiently accurate 
for the analytic purpose that motivates the inquiry. In other words, the simulation is 
terminated the first time a confidence interval width for a desired quantity is sufficiently 
small. We refer to such a procedure as a sequential fixed-width stopping rule and note the 
total simulation effort will be random. 

As we show later, fixed-width methods are especially desirable because they are theo- 
retically justified and constrained by few assumptions. The simplest fixed-width rule, first 



studied in MCMC by Jones et al. (2006), stops the simulation when the width of a confi- 



dence interval based on an ergodic average is less than a user-specified value, say e. Flegal 



et al. (2008) and Jones et al. (2006) show this stopping rule is superior to using convergence 



diagnostics as a stopping criteria. 

In this paper, we introduce relative fixed-width stopping rules that eliminate the need to 
specify an absolute value for e. Specifically, the simulation is terminated the first time the 
width of a confidence interval is sufficiently small relative to the size of a target parameter. 
We consider two measures of size, magnitude and standard deviation. Further, we illustrate 
the utility of relative fixed-width stopping rules for simultaneous estimation of multiple 
parameters. 

Specificity requires some notation. Let tt denote a probability distribution having sup- 
port X C R d , d > 1, about which we wish to make inference. This inference is usually based 
on some feature of tt denoted 9 n . For example, we may want a quantile of tt or if g : X — > K, 
we may need to calculate 



g(x)ir(dx) . 



Frequently tt is such that MCMC is the only viable technique for estimating 9 n . The 
basic MCMC method entails constructing a time-homogeneous Harris ergodic Markov chain 
on state space X with a- algebra B = B(X) and invariant distribution 
it. The popularity of MCMC methods result from the ease with which X can be simulated 
( Robert and CasellaJ |2004l ). 

Suppose we simulate X for n iterations, where n is finite. Define 9 n as an estimator of 
9 V from the observed chain. Outside of toy examples, no matter how long our simulation, 
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there will be an unknown Monte Carlo error, 9 n — 9 n . While it is impossible to assess this 
error directly, we can obtain its approximate sampling distribution if a Markov chain central 
limit theorem (CLT) holds. That is, if 



^[On-e^ A n(o,^) (i) 

as n — y oo where a^j G (0,oo). Denote Ag as the posterior variance associated with 9^. 
Then it is important to note that due to the correlation present in a Markov chain cr| 7^ Ag, 
except in trivial cases. 

For now, suppose we have an estimator such that a\ — > <t| almost surely as n — > 00. 
This allows construction of a (1 — 5)100% confidence interval for 6 n with width 

w s = 2z S / 2 ^ (2) 



where z$/2 is a critical value from a standard Normal distribution. The width at ([2]) al- 
lows analysts to report the uncertainty in their estimates and users to assess the practical 
reliability. Moreover, we will use ws to construct sequential fixed-width stopping rules. 

Our work advocates stopping the simulation the first time ws is sufficiently small. We 
consider three distinct stopping rules: (i) an absolute precision rule that terminates when 
wg < e, (ii) a relative magnitude rule that terminates when w$ < e\9^\ and (iii) a relative 
standard deviation rule that terminates when ws < eAg. 



The theoretical properties of (i) and (ii) have been studied by Glynn and Whitt (1992), 
which we extend to establish conditions for asymptotic validity of (iii). Asymptotic validity 
is important since it implies the simulation will terminate w.p.l and the resulting confidence 
intervals will have the right coverage probability. 



Flegal et al. (2008), Flegal and Jones (2010) and Jones et al. (2006) have previously 



investigated (i) for MCMC expectation estimation. We are not aware of any prior use of 
fixed-width methods for quantile estimation or any use of (ii) or (iii) as a stopping rule in 
MCMC. The rule (iii) has significant promise in Bayesian applications since the simulation 
terminates the first time the length of a confidence interval is less than an eth fraction of the 
magnitude of the standard deviation of 9 n . In other words, the simulation stops when an 
estimate of 9^ is sufficiently accurate relative to an associated posterior standard deviation. 
Another substantial benefit of rule (iii) is it easy to implement in multivariate settings since 
e can remain constant. 

There are two main assumptions for asymptotic validity. First, we require a limiting dis- 
tribution for the Monte Carlo error such as at ([!]). Second, we require a strongly consistent 
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estimator of the associated asymptotic variance, that is &\ — > a^j almost surely as n — > oo. 
We later discuss these assumptions in detail for estimating expectations and quantiles. 

Finally, we investigate the finite sample properties of relative fixed- width stopping rules 
through three examples. Our first example considers an independence Metropolis sampler to 
explore an exponential random variable. Our second example considers exploring a mixture 
of bivariate Normal distributions with Metropolis Hastings and Gibbs samplers. While 
these are only toy examples, we will use true parameter values to illustrate the utility of 
our stopping rules. Our final example considers a Bayesian version of a logistic regression 
to model the presence or absence of the freshwater eel Anguilla australis. 

Using these examples, we terminate the simulation with the three distinct fixed- width 
stopping rules and calculate confidence intervals for a vector of target parameters. Over 
replicated simulations, all the finite sample empirical coverage probabilities are close to a 
specified nominal level. Thus, fixed-width stopping rules provide a theoretically valid and 
practically accurate procedure to determine when to stop a MCMC simulation. 

For Bayesian practitioners, we advocate the relative standard deviation fixed-width 
stopping rule (iii) since it is easy to implement and applicable in multivariate settings 
without a priori knowledge of the target parameter size. As our examples show, setting 
e = 0.02 provides excellent results in a wide variety of univariate and multivariate settings. 

The rest of this paper is organized as follows. Section [2] formally introduces relative 
fixed- width stopping rules and establishes asymptotic validity. Section [3] investigates fixed- 
width stopping procedures when estimating expectations and quantiles. Section [4] studies 
the finite sample properties in three numerical examples and concludes with a discussion 
that provides some recommendations to practitioners. 

2 Sequential fixed-width procedures 

In this section, we obtain conditions that ensure asymptotic validity of fixed-width proce- 
dures. The primary assumption is a limiting distribution holds for the Monte Carlo error 
such as at Q. Actually, we require the limiting process at ^ satisfy a stronger condition, 
in particular, a functional central limit theorem (FCLT). 

Define an interval 



If ([I]) holds and a n is weakly consistent for erg, then C[n] achieves the nominal coverage 
level as the sample size n — > oo. Thus we have a valid confidence interval provided the 
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sample size is permitted to go to oo. 

Now consider a sequential procedure that terminates the simulation when the length of 
a confidence interval drops below a prescribed level e. We will refer to this type of stopping 
rule as an absolute precision fixed-width stopping rule. For such a rule, the time at which 
the simulation terminates is defined by 

T(e) = inf {n > : 2z$/ 2 &n/Vn < e} • 
Unfortunately, use of this stopping rule is insufficient because T(e) can terminate much too 



early if a n is poorly behaved for small n (Glynn and Whitt, 1992). Instead, suppose p(n) 
is a positive function that decreases monotonically such that p{n) = o(n -1 / 2 ) as n — > oo 
and let n* be the desired minimum simulation effort (a reasonable default is p(n) = el{n < 
n*) + Then an absolute precision stopping rule terminates the simulation at 



Ti(e) = inf [n > : 2z S / 2 v n /Vn + p( n ) < e} 



The following result, an immediate consequence of Theorem 1 in Glynn and Whitt 



(1992), yields asymptotic validity of the sequential stopping rule 7i(e). Note the desired 
coverage probability will be obtained in an asymptotic sense as e — > 0. 



Proposition 1. Suppose a FCLT at ([!]) holds. If ' a n — > oq w.p.l as n 

the simulation will terminate w.p.l and 



oo, then as n 



or e 



Pr G C[Ti(e)]) 



Remark 1. Glynn and Whitt (1992) show weak consistency of a n is not enough to ensure 



asymptotic validity . 

The stopping rule 71(e) has previously been used for estimating expectations in MCMC 



(Flegal et al. 2008 Flegal and Jones, 2010 Jones et al. 2006). We further show this rule 



works well for MCMC estimation of quantiles in the following section. The challenge in 
both settings is finding a strongly consistent estimator of o~q. 

One can consider a variant of the stopping rule T\{e) known as a relative precision 
stopping rule, which avoids having to choose an absolute value for e. Simply put, the 
simulation is run until the length of a confidence interval is less than an eth fraction of 
the magnitude of the parameter of interest, 9^. Using 6 n as an estimator of 9 n yields the 
following relative magnitude stopping rule 

T 2 (e) = inf |n > : 2z & i 2 o- n l ^ + p{n) < e 9 n j . 
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For large n, ?2(e) will behave like Ti(e\9 n \). The following obtains asymptotic validity of 



T2(e), which is a direct consequence of Theorem 3 in Glynn and Whitt (1992). 



Proposition 2. Suppose a FCLT at ([!]) holds and \9 T \ > 0. If n — > 9^ w.p.l and a n — > ag 
w.p.l as n —> oo, then as n — >■ oo or e — > the simulation will terminate w.p.l and 

Pr (9 n G C[T 2 (e)]) -> 1 - <5 . 

Note that Proposition [i] requires # ra — )• w.p.l along with necessary conditions of 
Proposition [TJ In general stochastic simulations, this condition does not immediately follow 
from ([!]) (see Example 2 of Glynn and Whitt , 1988 ) but is readily available when W is an 



expectation via the Markov chain strong law of large numbers (SLLN). 

While ?2(e) has some support in the operations research literature, it makes little in- 
tuitive sense in Bayesian settings. Specifically, if 9 V = then T2(e) will be theoretically 
invalid and poorly behaved in finite simulations. In addition, T^ie) could be problematic 
even when 9 n is close to zero, which we illustrate through example in Section [4} 

Given the popularity of MCMC in Bayesian settings, it is useful to consider another 
specifically designed variant of ?i(e). To this end, we propose a stopping rule that terminates 
the simulation when the length of a confidence interval is less than an eth fraction of the 
magnitude of Xg, i.e. the posterior standard deviation of 9 n . Suppose X n is an estimator of 
Xg and consider the following stopping rule 

T 3 (e) = inf |n > : 2z 5 / 2 a n /^n + p(n) < eA n | . 

For large n, T^{e) will behave like T\(eXg). The benefit of using Ts(e) is that e is selected 
as a fraction rather than in the units of the target parameter. Hence, a single value of 
e would be appropriate for target parameters of any magnitude. Naturally, decreasing e 
would decrease the uncertainty of the resulting estimates and could be done simultaneously 
for multiple parameters. The following establishes asymptotic validity of 23(e), which we 
prove in Appendix [A| 

Theorem 1. Suppose a FCLT at Q holds and Xg > 0. If X n — > Xg w.p.l and o n — > o~g 
w.p.l as n — >■ 00, then as n —> 00 or e — > the simulation will terminate w.p.l and 

Pr(6 v eC[T 3 (e)])-H-5. 

Note the only additional condition required for Theorem [l] is a strongly consistent esti- 
mator of Xg. For expectations, an estimator is readily available via the Markov chain SLLN. 
In the case of quantiles, we discuss a viable estimator in the following section. 
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The benefit of the stopping rule Ts(e) is twofold. First, one only needs to specify a rela- 
tive e, and hence no knowledge about the magnitude is required. Second, when estimating 
multiple parameters a single e will suffice to obtain estimates whose uncertainty will be 
comparable relative to their standard deviations. In other words, we have developed a sim- 
ple, yet informative, stopping criteria applicable in multivariate settings. In these settings, 
one could address the issue of multiplicity by adjusting the critical value appropriately. We 
illustrate this procedure via examples in Section [4j and show the resulting simultaneous 
confidence regions obtain at least the nominal coverage probability. 

Remark 2. Asymptotic validity of relative stopping rules can be established when ([!]) is 
replaced by a more general M- valued stochastic process (Glynn and Whitt, 1992[ ). The 
generalization enables consideration of 6 n that follow non-Normal asymptotic distributions. 



3 Applications 



This section demonstrates that fixed-width stopping rules are appropriate for MCMC es- 
timation of expectations and quantiles. This is an important contribution since we know 



of no other formal stopping criteria applicable in both settings. Raftery and Lewis (1992) 



propose a heuristic approach to terminating an MCMC simulation when the primary inter- 



est is quantile estimation. However, Brooks and Roberts (1999) argue "in the case where 



quantiles themselves are not of interest, this method should be used with caution". 

First, we require a bit more notation to describe sufficient mixing conditions for a Markov 
chain CLT and consistent estimation of the asymptotic variance. An interested reader is 



directed to Meyn et al. (2009) and Roberts and Rosenthal (2004) for more on Markov chain 
theory. 

Recall A is a Harris ergodic Markov chain on state space X with u-algebra B = B(X) and 
invariant distribution ir. Denote the n-step Markov kernel associated with A as P n (x,dy) 
for n G N. Then if A € B(X) and k G {0, 1, 2, . . .}, P n (x, A) = Pr(A fc+n G A\X k = x). Let 
|| • || denote the total variation norm. Let M : X >— > M + and 7 : N 1— > M + be decreasing such 
that 

\\P n (x r )-ir(.)\\<M(x) 7 (n). (3) 

Polynomial ergodicity of order m where m > means ((3j) holds with E W M < 00 and 
7(71) = n~ m for all Ao = x. Geometrical ergodicity means (|3j) holds with 7(71) = t n for 
some < t < 1 for all Xq = x. Uniform ergodicity means ^ holds with M bounded and 
7(71) = t n for some < t < 1. 

Establishing ^ directly can be challenging, but some constructive techniques are avail- 



able (Jarner and Roberts 2002 Meyn et al. 2009). Most literature on MCMC algorithms 
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focuses on establishing geometric and uniform ergodicity, see e.g. Hobert (2011), Jones and 



Hobert (2001), Johnson et al. (2011), Mengersen and Tweedie (1996), Roberts and Tweedie 



(1996) and Tierney (1994). Less has been said concerning polynomial ergodicity, but an 



interested reader is directed to Douc et al. (2004), Fort and Moulines (2000), Fort and 



Moulines (2003), Jarner and Roberts (2002), Jarner and Roberts (2007) and Jarner and 



Tweedie (2003). 



3.1 Expectations 



For MCMC estimation of an expectation, one can obtain all the necessary conditions for 
asymptotic validity of fixed-width stopping rules. Let g : X — > K, then we consider estima- 
tion of 



fig := E n [g(X)] = / g(x)n(dx) . 
Jx 

Estimating \i g is natural by appealing a Markov chain SLLN, a special case of the Birkhoff 



Ergodic Theorem (p. 558 Fristedt and Gray, 1997). Specifically, if E n \g\ < oo then w.p.l 



9n 



. n— 1 



n 



8=0 



g{x 



(ih 



(jL g as n — > oo 



Hence the SLLN yields strongly consistent estimators of fx g and X 2 , = Var[g] (provided 
Eng 2 < oo) necessary for Proposition [2] and Theorem [TJ respectively. 

We can obtain an approximate sampling distribution for the Monte Carlo error via a 
Markov chain CLT if 



V^(5n-^ g )4 N(0,a 9 2 ) 



(4) 



as n —7- 00 where a g G (0, 00). Conditions that ensure Q can be found in Chan and Geyer 
( |1994[ ), |Jonesl ( |2004D , [Meyn et aI.| ( |2009| |, |Roberts and Rosenthal| fl2004[ ) and |Tierney| ( |1994[ ). 
For example, if X is geometrically ergodic and E 1T \g\ 2+e < 00 for some e > 0, then Q holds. 
Fortunately, Markov chains frequently enjoy a FCLT under the same conditions ( jlbragimov 



1962 Oodaira and Yoshihara, 1972) 



There are many strongly consistent variance estimation techniques applicable for a g 



in MCMC settings including batch means (Flegal and Jones 2010 Jones et al. 2006), 



spectral variance techniques (Flegal and Jones, 2010) and regenerative simulation (Hobert 



et al. 2002, Mykland et al. 1995). We consider only non-overlapping batch means (BM) 



because it is easy to implement and available in many software packages, e.g. the mcmcse 
package available on CRAN. 

In BM the output is broken into a n batches where each batch is b n iterations in length. 
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Suppose the algorithm is run for a total of fi — d n b n iterations ctnci define 



Yi 



The BM estimate of is 



=C?-l)6n+l 



for j = 1, 



b n 



(5) 



J'=l 



In general, the BM estimator at §5§ is not a consistent estimator of cr?. However, 



Jones 



et al. 



(2006) establish necessary conditions for <r„ 



0"^ with probability 1 as n 



oo if 



the batch size and the number of batches are allowed to increase as the overall length of 
the simulation increases. Setting b n = [n T \ and a n = [n/b n \, the regularity conditions 
require that X be geometrically ergodic, E 7r \g\ 2+ei+£2 < oo for some e% > 0, €2 > and 
(1 + ei/2) _1 < r < 1. A common choice of r = 1/2 (i.e., b n = [\/n\ and a n = [n/b n \) has 



been shown to work well in applications (Flegal et al. 2008 |Flegal and Jones 2010, Jones 



et al. 2006). 



Remark 3. Most sampling plans require storing the entire Markov chain to allow for re- 
calculations as the batch size increases with n. If storage is a concern, one could consider 
increasing the batch size of the form b n £ {2, 4, 8, 2 k , ...} in an effort to reduce memory 
usage. One can establish strong consistency for the BM variance estimator with such a 



sampling plan using results in Jones et al. ( 2006 ) and Bednorz and Latuszynski ( 2007 ) . 



3.2 Quantiles 

It is routine to estimate univariate quantiles associated with tt, especially in Bayesian ap- 
plications. To this end, let W ~ tt and recall g : X — >• M. Setting V = g(W), we consider 
estimation of the quantiles associated with the univariate distribution of V. Suppose Fy 
denotes the cumulative distribution function of V, then our goal is to obtain 

i q ■= Fy\q) = ku> : F v (v) > q} . 
Little has been formally said regarding MCMC estimation of quantiles, but we outline the 



current state of understanding (for more details see Flegal et al. 2012). 

A natural estimator of £ q is the inverse of the empirical distribution function given by 

L,q ■= Y n(j+i) where 3 <nq<j + 1 , (6) 
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where Y n fj\ denotes the jth order statistic of {Yq, . . . , Y n -\} = {g(Xo), . . . , g(X n _i), }. If X 
is Harris recurrent and £ g is the unique solution y of Fy{y—) < q < Fy(y), then £ n)(J — > £ q 
w.p.l as n — > oo (Flegal et al. 2012). 

Under stronger mixing conditions on X, one can obtain a Markov chain CLT. To this 
end, define 

oo 

a 2 {y) := Var, [I(Y < y)\ + 2 £ Cov n [I(Y < y), I(Y k < y)\ . 

k=l 

Suppose there is e > such that X is polynomially ergodic of order 2.5 + e. If Fy has a 
density fy positive and bounded in some neighborhood of £„, then as n — > oo 



(7) 



where 7 2 (£ ? ) = a 2 (t q )/[fv{Z q )} 2 , provided a 2 (£ q ) > ( |Flegal et al.| |2012[ ). A FCLT is 
extremely likely to hold under similar conditions, but we are unaware of a formal proof. 

Estimation of the variance from the asymptotic Normal distribution at ([7]) is broken 
into two parts. First, we plug in £ nj(3 for £ q and separately consider estimating fy{^ n ^ q ) and 
c" 2 (Cn,(j)- Estimating fy{^ n ,q) uses a kernel density approach with a gaussian kernel, which 
we denote as fv(£n,q)- There are well known conditions guaranteeing strongly consistent 
estimation of the density at a point (see e.g. Kim and Lee, 2005, Yu, 1993). 

We will use BM for estimating <r 2 (^ n . g ). Suppose we have n = a n b n iterations, then 
for k = 0, . . . ,a n - 1 define Uk(€n,q) ■= K l H^q 1 I( Y kb n +i < £n,q)- The BM estimate of 



^ a-n-1 

^BM(in,q) = — 7 (Uk(in,q) ~ £Ai(£n,< 

a n J- V 
fe=0 

Combining fy(U,q) and <r| M (Cn,q), we estimate 7 2 (£ g ) with 



BM 



(£n )9 ) 



fv(Cn,q) 

This approach is implemented in the R package mcmcse which is used to perform the com- 



putations in our examples. Flegal et al. (2012) outline the conditions that ensure strong 
consistency of this estimator. 

The relative standard deviation fixed- width stopping rule of Theorem [T] requires a esti- 

q(l-q) 



mation of 



Xe 



Mtq 
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We use the same kernel density estimate resulting in 



g(l - q) 

fv(€n,q) 



4 Numerical studies 



This section investigates the finite sample properties of fixed- width stopping rules through a 
variety of simulations. In each example, we independently repeat the MCMC simulation to 
evaluate the resulting finite sample confidence intervals. Naturally, this evaluation requires 
the true parameter values. In our first two examples, the true values are readily available. 
In our final example, the truth was estimated using an independent long run of the MCMC 
sampler. Overall, the empirical coverage probabilities obtained via fixed-width stopping 
rules are remarkably close to the nominal level. 

Each simulation considered both expectations and quantiles with the following common 
methodology. For a single replication, the same MCMC draws were used in applying the 
three stopping rules. Further, we uniformly set p(n) = el(n < n*) + and estimate an 
via BM methods with b n = [\/n\ calculated with the mcmcse package. Finally, standard 
errors for the empirical coverage probabilities equal y/p(l — p)/r where r is the number of 
replications. 



4.1 Exponential distribution 

Consider an Exp(l) target distribution, i.e. f(x) = e~ x I(x > 0). It is easy to show that 
E[X] = 1 and F~ 1 (q) = log(l - g) _1 , which we use to evaluate finite sample confidence in- 
tervals obtained via fixed- width methods. We will sample from f(x) using an independence 
Metropolis sampler with an Exp(l/2) proposal and note this chain is geometrically ergodic 



(Jones and Hobert, 2001). 



First, consider estimation of E[X] using each combination of Tj(e) for i E {1,2,3} and 
e G {0.10,0.05,0.02}. The chain was started from 1 and ran for a minimum of n* = 1000 
iterations. If the stopping criteria was not met, an additional 500 iterations were added 
to the chain before checking again. The simulation was repeated for 2000 replications to 
evaluate the resulting coverage probabilities. 

Table [T] summarizes the mean and standard deviation of the number of iterations at 
termination along with the resulting coverage probabilities. All the coverage probabilities 
are close to the 0.90 nominal level suggesting all three stopping rules are preforming well. 
Note the mean iterations are approximately equal, which is expected since E[X] = 1 and 
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\ e = VarpT] = 1. 





Length (SD) E[X] 


Length (SD) £ 5 


7i(0.10) 
71(0.05) 
Ti(0.02) 


2.44E3 (4.9E2) 0.884 
8.89E3 (1.2E3) 0.894 
5.36E4 (4.7E3) 0.887 


2.70E3 (5.9E2) 0.858 
1.01E4 (1.5E3) 0.881 
6.17E4 (5.4E3) 0.877 


r 2 (o.io) 

T 2 (0.05) 
T 2 (0.02) 


2.44E3 (4.8E2) 0.889 
8.90E3 (1.2E3) 0.891 
5.35E4 (4.7E3) 0.887 


5.40E3 (9.4E2) 0.880 
2.07E4 (2.4E3) 0.882 
1.29E5 (9.1E3) 0.883 


r 3 (o.io) 

T 3 (0.05) 
T 3 (0.02) 


2.45E3 (4.7E2) 0.888 
8.90E3 (1.2E3) 0.888 
5.35E4 (4.6E3) 0.889 


2.79E3 (5.2E2) 0.865 
1.03E4 (1.3E3) 0.882 
6.23E4 (5.2E3) 0.877 



Table 1: Summary of coverage probabilities for estimation of E[X] and £.5 based on 2000 
replications and 0.90 nominal level. 

Next, consider estimation of the median, £5, using the same simulation settings. Table[l] 
summarizes the results from 2000 replications. Again the results are very close to the 0.90 
nominal level, though slightly lower than those for estimating the mean. Here we have 
£.5 = 0.693 and 0.5(1 — 0.5)/e^- 5 = 1, hence for fixed e we expect T\(e) and 13(e) to be 
similar and T 2 (e) to be larger. 

Finally, consider estimating the mean and an 80% Bayesian credible region simultane- 
ously, which we denote as $ = (E[X],£_i,£$). Due to increased computation time, each 
chain was run for a minimum of n* = 10000 iterations with an additional 5000 added 
between checks. The simulation was terminated the first time the length of a confidence 
interval was sufficiently small for each parameter in $. To adjust for multiplicity, we apply a 
Bonferonni approach. Specifically, we set individual confidence intervals to have a coverage 
probability of 0.90 1 / 3 = 0.9655 resulting in a simultaneous confidence region with coverage 
probability of at least 0.90. 

The simulation was repeated for 2000 replications with each combination of Tj(e) for 
i G {1,2,3} and e G {0.10,0.05,0.02}. Tabled summarizes the simulation results. We can 
see the individual coverage probabilities improve as e decreases, especially in the case of £.1. 
For e = 0.02, all the individual coverage probabilities are remarkably close to the nominal 
level of 0.9655. Note the observed confidence region coverage probabilities are above the 
0.90 nominal level, which is unsurprising due to correlation between parameters in <3?. 
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Length (SD) 


E\X\ 


F i 
s.l 


f n 


T? pci nn 

J_ IL^lUll 


±l{y.iv) 


Zr.oorj'4: ^o.yrjj j 


u.yuo 




U.c/UO 




J- 1 ^U.Ud J 


1 f)7F^ /'Q 7F3^ 




fl Q7Q 






Ti (0.02) 


6.53E5 (3.3E4) 


0.965 


0.967 


0.968 


0.917 


T 2 (0.10) 


6.71E4 (5.9E3) 


0.969 


0.979 


0.964 


0.925 


T 2 (0.05) 


2.29E5 (1.4E4) 


0.966 


0.974 


0.963 


0.920 


T 2 (0.02) 


1.29E6 (5.0E4) 


0.964 


0.963 


0.970 


0.915 


r 3 (o.io) 


1.00E4 (0) 


0.962 


0.991 


0.955 


0.927 


T 3 (0.05) 


2.31E4 (2.9E3) 


0.963 


0.983 


0.958 


0.921 


T 3 (0.02) 


1.30E5 (9.1E3) 


0.961 


0.970 


0.965 


0.914 



Table 2: Summary of coverage probabilities for estimation of based on 2000 replications. 
Individual confidence intervals have a 0.9655 nominal level, resulting in a 0.90 nominal level 
confidence region. 



4.2 Mixture of bivariate Normals 

Consider a mixture of bivariate Normals X 



X u X 2 f = pYi + (1 - p)Y 2 , where 



Vu" 


~N 2 ( 


Mil 


3 


"*li o " 




Y 12 _ 




A*12 




a\ 2 _ 


J and 



v 2 r 


~N 2 ( 


M21 




*21 





y 22 _ 




M22 










In this example, we choose p = 0.25, /in = 1, [i\ 2 = 10, /i 2 i = 2.5, \x 22 = 25, an = 0.5, 
cri2 = 5, <7 2 i = 0.7 and cr 22 = 7. 

We first sample from /(X) with two different component- wise Metropolis random walk 
algorithms, one with Uniform proposals and another with Normal proposals. For the Uni- 
form proposals, we apply a Unif(— 3,3) and Unif(— 30, 30) random walk for the X\ and 
X 2 dimensions, respectively. For the Normal proposals, we apply a iV(0,3 2 ) and iV(0,30 2 ) 
random walk for the X\ and X 2 dimensions, respectively. It can be shown that these chains 
are geometrically ergodic ( Jarner and Hansen , |2000 ) . 

Consider estimation of = (E[X],£.i,£$) using fixed-width stopping rules Tj(e) for 
i € {1,2,3} and e G {0.10,0.05,0.02}. We ran the chain for a minimum of n* = 5000 
iterations and added 1000 iterations between checking the stopping criteria. This simulation 
was repeated for 1000 independent replications. 

Table [3] summarizes the mean and standard deviation of the number of iterations at 
termination along empirical coverage probabilities from the Uniform and Normal proposals. 
Notice for both samplers, the coverage probabilities improve as e decreases and are close to 
the 0.95 nominal level once e = 0.02. It appears the Metropolis random walk with Normal 
proposals is mixing faster since the overall simulation effort is substantially lower than that 
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of the Uniform proposals. This difference in simulation effort illustrates the importance of 
specifying a good proposal distribution in MCMC simulations. 



Uniform 


Length (SD) 


E[Xi] C.9,Xi E[X 2 ] £.i,x 2 £,.9,x 2 


Ti(0.10) 
Ti (0 05) 

J- I l KJ * KJ KJ 1 

71(0.02) 


14,658 (3.4E3) 
59 869 (9 1E3) 

KJ kJ y KJ \J kJ \ \J • -M- -■ — J KJ J 

391,566 (3.1E4) 


0.930 0.932 0.917 0.936 0.945 0.937 
934 922 939 940 934 953 

KJ • \J KJ -I- KJ • kJ £J £J \J • KJ kJ \J • %J J- V_/ W • %J KJ -l- W • y ".7 

0.956 0.944 0.945 0.956 0.948 0.953 


r 2 (o.io) 

Ti(0 05) 
T 2 (0.02) 


20,897 (5.0E3) 
85 401 (1 2E4) 

KJ KJ ^ -L V ' J_ \ -L ■ * — J -i- J 

556,821(3.9E4) 


0.929 0.933 0.911 0.931 0.936 0.938 
950 926 934 929 925 942 
0.953 0.946 0.954 0.950 0.938 0.956 


T 3 (0.10) 
Ti(0 05) 
T 3 (0.02) 

Normal 


8,827 (1.0E3) 
35 733 (2 9E3) 

KJ KJ ^ 1 KJ KJ \ ^ • kJ -I — i KJ f 

233,312 (1.3E4) 
Length (SD) 


0.926 0.928 0.899 0.920 0.922 0.920 
924 938 931 934 928 937 

V / • \J — -i- V / • ' / KJ KJ V J • ■ KJ J- V t • \ ' KJ -l V / • ■ / £J KJ KJ • \J KJ 1 

0.954 0.955 0.959 0.948 0.958 0.956 

E _X 1 C1 ^ C Q -A^9 CI C Q 

L J- J s.i,Ai s.y,vvi i £\ s-J-j-^-2 s.y,^v2 


Ti(0.10) 
71(0.05) 
71(0.02) 


8,028 (1.5E3) 
29,844 (3.7E3) 
186,061 (1.3E4) 


0.946 0.939 0.939 0.934 0.943 0.937 
0.927 0.936 0.948 0.917 0.932 0.953 
0.952 0.936 0.952 0.943 0.946 0.938 


T 2 (0.10) 
T 2 (0.05) 
T 2 (0.02) 


11,307 (2.1E3) 
42,338 (4.6E3) 
261,741 (1.6E4) 


0.949 0.933 0.940 0.940 0.944 0.943 
0.911 0.943 0.956 0.937 0.934 0.951 
0.940 0.938 0.956 0.949 0.938 0.945 


T 3 (0.10) 
T 3 (0.05) 
T 3 (0.02) 


5,114 (3.2E2) 
17,654 (1.8E3) 
112,626 (7.6E3) 


0.944 0.950 0.933 0.936 0.936 0.924 
0.922 0.930 0.943 0.925 0.921 0.939 
0.933 0.946 0.941 0.941 0.930 0.940 



Table 3: Summary of coverage probabilities for estimations of <3? using a Metropolis random 
walk with Uniform and Normal proposals based on 1000 replications and a 0.95 nominal 
level. 



Next, we consider a Gibbs sampler using the full conditional densities, i.e. 

/x^^iM = Px 2 Y n + (1 - Px 2 )Y 2 i and 
fx 2 \ Xl (x2\xi) = P Xl Y 12 + (1 - P Xl )Y 22 , 



where 



(l-p)ai2 j 1 / (x 2 - Aii2\ 2 (x 2 -^ 22 x2 



p x 2 = H exp 



P&22 \ ^ \ \ Cl2 / \ <J 22 



and 



n ( -i i (1-P)0"11 ] 1 / /Xl — /Xn\ (xi-^21 

Px 1 = H exp 



PO"21 M \ V °"H / V °"21 
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Note, X1IX2 = %2 and X 2 |Xi = x\ are easy to sample from since they are mixtures of 
Normal random variables. 

Table [4] summarizes the results for the Gibbs sampler. Notice, the coverage probabilities 
do not uniformly improve as e decreases. However, they are all close to the nominal 0.95 
level using significantly fewer total iterations, suggesting the Gibbs sampler mixes better 
than either of the Metropolis random walk samplers. 

As a final comparison, we performed additional simulations via i.i.d. sampling (not 
shown). The resulting empirical coverage probabilities were similar to using the Gibbs 
sampler, albeit with slightly fewer iterations. 



Gibbs 


Length (SD) 


£LYi] 




£.9,Xi 


E[X 2 ] 


tl,X 2 


£,.9,X 2 


Ti(0.10) 


1,930 (3.7E2) 


0.941 


0.940 


0.937 


0.954 


0.958 


0.927 


21(0.05) 


5,727 (8.7E2) 


0.946 


0.958 


0.941 


0.942 


0.945 


0.940 


71(0.02) 


31,170 (2.8E3) 


0.935 


0.945 


0.961 


0.937 


0.937 


0.944 


r 2 (o.io) 


2,465 (5.4E2) 


0.935 


0.939 


0.939 


0.954 


0.950 


0.937 


T 2 (0.05) 


7,865 (1.1E3) 


0.950 


0.959 


0.943 


0.955 


0.954 


0.952 


T 2 (0.02) 


43,756 (3.6E3) 


0.933 


0.936 


0.959 


0.936 


0.959 


0.946 


T 3 (0.10) 


1,182 (3.9E2) 


0.929 


0.936 


0.942 


0.936 


0.936 


0.924 


T 3 (0.05) 


3,786 (6.2E2) 


0.956 


0.951 


0.944 


0.940 


0.940 


0.935 


T 3 (0.02) 


20,289 (2.0E3) 


0.945 


0.947 


0.954 


0.940 


0.943 


0.952 



Table 4: Summary of coverage probabilities for estimations of using a Gibbs sampler 
based on 1000 replications and a 0.95 nominal level. 



4.3 Bayesian logistic regression 



Our final example considers the Anguilla eel data provided in the dismo R package (see e.g. 



Elith et al. 2008 Hijmans et al. 2010). The data consists of 1,000 observations from a New 



Zealand survey of site-level presence or absence for the short-finned eel (Anguilla australis). 
We selected six out of twelve covariates as in |Leathwick et al.| ( |2008| ). Five are continuous 
variables: SegSumT, DSDist, USNative, DSMaxSlope and DSSlope; one is a categorical 
variable: Method, with five levels Electric, Spo, Trap, Net and Mixture. 

Let Xi be the regression vector of covariates for the ith observation of length k and 
ii = (j3o, . . . ,j3g) be the vector regression coefficients. For the ith observation, suppose 
Yi = 1 denotes presence and Yi = denotes absence of Anguilla australis. Then the 
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Bayesian logistic regression model is given by 



Yi ~ Bernoulli{jpi) 

exp(xf/3) 
Pi ~ 1 + exp(xf 0) 

^~iV(0,o-|l fc ) , 



and, 



where is the x k identity matrix. For the analysis, cr| = 100 was chosen to represent a 
diffuse prior distribution on /? (Boone et al. , 2012). Further, we use the MCMClogit function 



in the MCMCpack package to sample from the target Markov chain. 

Suppose we are interested in estimating the posterior mean along with an 80% Bayesian 
credible interval for each regression coefficient in the model. Given that we are working with 
real data, the true values are naturally unknown. Instead, we ran 1000 independent chains 
for 1E6 iterations to obtain an accurate estimate, which we treat as the truth (Table pjl). 



Variable 


Pi 


& 


Ai) 

S.9 


Intercept 


-10.463 (2.7E-5) 


-12.224 (3.9E-4) 


-8.730 (3.7E-4) 


SegSumT 


0.657 (1.5E-5) 


0.559 (2.1E-5) 


0.757 (2.2E-5) 


DSDist 


-4.02E-3 (3.3E-7) 


-6.15E-3 (4.9E-7) 


-1.93E-3 (4.4E-7) 


USNative 


-1.170 (7.1E-5) 


-1.625 (9.9E-5) 


-0.718 (1.0E-4) 


MethodMixture 


-0.468 (6.8E-5) 


-0.910 (9.8E-5) 


-0.028 (9.8E-5) 


MethodNet 


-1.525 (8.2E-5) 


-2.026 (1.2E-4) 


-1.035 (1.1E-4) 


MethodSpo 


-1.831 (1.3E-4) 


-2.623 (2.2E-4) 


-1.798 (1.4E-4) 


MethodTrap 


-2.594 (1.1E-4) 


-3.285 (1.8E-4) 


-1.937 (1.3E-4) 


DSMaxSlope 


-0.170 (1.1E-5) 


-0.244 (1.7E-5) 


-0.099 (1.5E-5) 


USSlope 


-0.052 (3.7E-6) 


-0.076 (5.5E-6) 


-0.028 (5.1E-6) 



Table 5: Summary of estimated true values with standard errors for the Bayesian logistic 
regression example. 



Consider estimating <3?j = f/3j,£.i j£.9 J for j = 0,...,9 using fixed-width stopping 
rules Tj(e) for i G {1,2,3}. From the magnitudes in Table [5j it is easy to see a single e 
would be problematic for Ti(e). Instead, we will specify an e for each $j with respect to its 
magnitude. Specifically, we choose three simulation settings such that €\ = (1, 0.01, 0.001, 
0.1, 0.1, 0.1, 0.1, 0.1, 0.01, 0.01), 0.5ei and 0.2ei. 

A single e value for T2(e) will also be problematic since there are parameters with very 
small absolute values (e.g. DSDist). We instead specify an e for each <F,-. In this case, we 
choose three simulation settings such that 62 = (0.1, 0.1, 1, 0.1, 1, 0.1, 0.1, 0.1, 0.1, 1), 0.5e2 
and 0.2C2- 
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For both Ji(e) and 22(e), it becomes overwhelmingly tedious to specify appropriate e 
vectors when the number of parameters becomes large. However, for the stopping rule 13(e) 
we can use a single e for the 30 dimensional target parameter vector. Specifically, we choose 
three simulation settings such that 63 E {0.10,0.05,0.02}. 

For the two larger e settings, we set n* = 10000 and added 1000 iterations between 
checks. For the smallest e setting, we set n* = 1E5 and added 10000 iterations between 
checks due to increased computational demands. Each simulation setting was repeated 1000 
times independently. 

Table [6] summarizes the empirical coverage probabilities. We can see the coverage prob- 
abilities for each stopping rule increase towards the nominal level of 0.95 as e decreases, 
suggesting that all the stopping rules perform well. For high dimensional settings such as 
this, 13(e) provides a distinct practical advantage since a practitioner can specify a single e 
value. 

To adjust for multiplicity, we again apply a Bonferonni approach. We set individual 
confidence intervals to have a nominal level of 0.80 1//10 = 0.9779 resulting in simultaneous 
confidence region with nominal level of at least 0.80. We only considered estimating the 
posterior mean of the 10 dimensional vector using T^(e) with e G {0.20,0.10,0.05,0.02}. 
The minimum simulation effort was n* = 1E5 iterations with an additional 1000 added 
between checks. Again, for the smallest e setting, we set n* = 1E6 with an additional 
10000 added between checks. The simulation was terminated the first time T^(e) was met 
and repeated 1000 times independently. 

Table [7] summarizes the simulation results. We can see that, as e decreases, all the 
individual coverage probabilities are remarkably close to the nominal level of 0.9779. Note 
the observed confidence region coverage probabilities approach the nominal level of 0.80 as 
expected. However, it is bit surprising how close this is to the nominal 0.80 level given 
possible correlation among parameters. To this end, we investigated the correlation be- 
tween pairs of target parameters. We found that most pairs have low correlation, except for 
strong correlation between (Intercept, SegSumT) and moderate correlation between (US- 
Native, USSlope). Given the lack of correlation, the confidence region coverages are very 
encouraging. 

4.4 Discussion 

This paper considers absolute precision, relative magnitude, and relative standard deviation 
fixed- width stopping rules in the context of MCMC simulations. Under limited assumptions, 
we show fixed-width stopping rules obtain a desired coverage probability in an asymptotic 
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ri(6i) 


Ti(0.5ei) 


Ti(0.2ei) 


Variable 


Pj 


Aj) 

S.i 


Aj) 

S.9 


Pj 


A'j) 
s.i 


Aj) 

S.9 


Pj 


Aj) 
s.i 


Ai) 

S.9 


Intercept 


0.936 


0.933 





912 


0.937 


0.942 





942 


0.946 


0.946 





930 


SegSumT 


0.932 


0.922 





916 


0.942 


0.941 





934 


0.953 


0.944 





936 


UbDlSt 


11.987 


0.969 





979 


11.975 


U.969 





960 


U.9oo 


U.954 





952 


yj oiy dLi ve 




0.929 





917 




0.933 





943 




0.939 





944 


MethodMixture 


0.930 


0.928 





920 


0.946 


0.948 





938 


0.935 


0.953 





940 


MethodNet 


0.946 


0.922 





936 


0.941 


0.948 





932 


0.943 


0.939 





935 


MethodSpo 


0.913 


0.913 





927 


0.931 


0.929 





931 


0.943 


0.942 





926 


MethodTrap 


0.928 


0.906 





937 


0.938 


0.930 





927 


0.941 


0.947 





947 


DSMaxSlope 


0.932 


0.930 





921 


0.942 


0.943 





945 


0.953 


0.958 





951 


USSlope 


0.921 


0.928 





935 


0.951 


0.927 





954 


0.957 


0.952 





962 


Length (SD) 


19, 


521 (3.8E3) 


76 


894 (9.5E3) 


492 


,910 (3.4E4) 





T 2 (e 2 ) 


T 2 (0.5e 2 ) 


T 2 (0.2e 2 ) 


Variable 


Pi 


Ai) 
s.i 


Aj) 

S.9 


Pi 


Ai) 
s.i 


Aj) 

S.9 


Pi 


Ai) 
s.i 


Aj) 

S.9 


Intercept 


0.928 


0.938 


0.915 


0.950 


0.948 


0.947 


0.945 


0.949 


0.938 


oegouml 


0.923 


0.916 


0.937 


0.953 


0.955 


0.948 


0.944 


0.947 


0.947 


USlJlSt 


0.985 


0.968 


0.975 


0.970 


0.958 


0.958 


0.956 


0.955 


0.947 


TTCJlVntiVp 
U OIN dLl ve 


0.921 


0.936 


0.921 


0.946 


0.933 


0.945 


0.940 


0.956 


0.941 


MethodMixture 


0.941 


0.938 


0.933 


0.942 


0.945 


0.916 


0.935 


0.933 


0.942 


MethodNet 


0.942 


0.920 


0.922 


0.940 


0.942 


0.939 


0.942 


0.944 


0.935 


MethodSpo 


0.919 


0.901 


0.924 


0.936 


0.923 


0.937 


0.947 


0.956 


0.947 


MethodTrap 


0.935 


0.910 


0.936 


0.939 


0.939 


0.931 


0.941 


0.933 


0.941 


DSMaxSlope 


0.937 


0.942 


0.916 


0.948 


0.942 


0.950 


0.942 


0.954 


0.955 


USSlope 


0.935 


0.933 


0.930 


0.949 


0.936 


0.941 


0.949 


0.944 


0.943 


Length (SD) 


37,667 (3.5E4) 


151,276 (8.9E4) 


1,161,400 (2.6E5) 






T 3 (0.10) 






T 3 (0.05) 






T 3 (0.02) 




Variable 


Pi 


& 


M 

S.9 


Pi 




Aj) 

S.9 


Pi 


& 


Aj) 

S.9 


Intercept 


0.932 


0.944 


0.929 


0.943 


0.950 


0.943 


0.943 


0.954 


0.934 


SegSumT 


0.932 


0.935 


0.941 


0.942 


0.934 


0.946 


0.942 


0.934 


0.946 


DSDist 


0.981 


0.969 


0.969 


0.968 


0.966 


0.955 


0.957 


0.954 


0.950 


USNative 


0.939 


0.942 


0.923 


0.941 


0.948 


0.954 


0.942 


0.943 


0.940 


MethodMixture 


0.939 


0.928 


0.920 


0.947 


0.943 


0.933 


0.927 


0.947 


0.928 


MethodNet 


0.929 


0.922 


0.931 


0.939 


0.939 


0.934 


0.930 


0.938 


0.939 


MethodSpo 


0.915 


0.902 


0.925 


0.924 


0.933 


0.926 


0.948 


0.946 


0.935 


MethodTrap 


0.930 


0.909 


0.920 


0.941 


0.937 


0.933 


0.939 


0.935 


0.948 


DSMaxSlope 


0.941 


0.932 


0.930 


0.940 


0.950 


0.943 


0.958 


0.955 


0.951 


USSlope 


0.939 


0.928 


0.940 


0.953 


0.937 


0.955 


0.954 


0.957 


0.958 


Length (SD) 


24,404 (1.4E3) 


78,886 (4.2E3) 


439,260 (1.7E4) 



Table 6: Summary of coverage probabilities for Bayesian logistic regression example with 
1000 independent replicates. The coverage probabilities have a 0.95 nominal level. 
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13{V.ZI)) 


l3(l).W) 


i3(U.U0J 




Variable 


a 

Pi 


a 

Pi 


a 

Pi 


a 

Pi 


Intercept 


u.yoy 


u.y t d 


n Q7fi 
u.y to 


u.y i o 


Qnn-QnmT 

oegoumi 


u.you 


u.y i i 


n Q7Q 

u.y < y 


u.y 1 4 


IJolJISTi 


u.yyo 


u.yoy 


u.yyo 


n Q7Q 

u.y < y 


u din aiive 


n Q/LR 


n Q78 
u.y i o 


D Q7D 

u.y / u 


n Q7^ 
u.y * o 


MpthnHMiYtnrp 


n 950 


973 


Q67 


Q68 


MethodNet 


0.962 


0.962 


0.976 


0.973 


MethodSpo 


0.946 


0.954 


0.968 


0.979 


MethodTrap 


0.950 


0.960 


0.970 


0.978 


DSMaxSlope 


0.966 


0.971 


0.977 


0.974 


USSlope 


0.964 


0.965 


0.973 


0.982 


Region 


0.693 


0.763 


0.792 


0.805 


Length (SD) 


10,082(2.7E2) 


29,729(1.8E3) 


100,261(5.2E3) 


583,488(1.9E4) 



Table 7: Summary of coverage probabilities for based on Ts(e) with 1,000 replicates. 
The coverage probabilities have a 0.9779 nominal level, resulting in a 0.80 nominal level 
confidence region. 

sense as e tends to 0. Moreover, we illustrate these rules perform well in a variety of finite 
sample settings provided e is specified to be small enough. 

A practical MCMC stopping rule should be applicable for a large number of parameters 
since practitioners usually report multiple expectation and quantile estimates. Unfortu- 
nately, choosing a single e could be problematic for absolute precision and relative magni- 
tude stopping rules. These stopping rules would be better served by specifying an e vector, 
which can be tedious when the number of parameters becomes large. 

Instead, we advocate use of the relative standard deviation stopping rule since it is 
easy to implement and applicable in multivariate settings without a priori knowledge of 
the target parameter size. Simply put, this rule terminates an MCMC simulation when 
estimates of target parameters are sufficiently accurate relative to their associated posterior 
standard deviations. The resulting estimates are approximately e _1 more accurate than 
their posterior standard deviations. We recommend using e = 0.02, which provided excellent 
results in the wide variety of examples considered here. However, a smaller e may be 
appropriate when the accuracy of estimation is critical. 

In any MCMC simulation, a key component is choosing a Markov chain that mixes 
well while sufficiently exploring the state space. As in the mixture of bivariate Normals, 
the sampler choice affects the performance significantly in terms of coverage probabilities. 
Moreover, the computational effort to achieve a reasonable accuracy varies depending on 
the sampling scheme. In practice, the true parameters values are unknown and thus poorly 
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behaved samplers may lead to suspicious inference. We have offered limited guidance in 
this direction, but note this is usually the most challenging aspect of an MCMC simulation. 



An interested reader is directed to Brooks et al. (2010) and the references therein for advice 
on sampling schemes. 

Finally, our examples only consider BM to estimate the asymptotic variance from a 
CLT since it is the most popular technique and widely available. Improving the variance 
estimation step might be possible using alternative methods such as overlapping batch 



means, spectral variance, or subsampling bootstrap methods (Flegal, 2012 Flegal and Jones 



2010 Flegal et al. 2012), which are currently available in the mcmcse package. 
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A Proof of Theorem [T] 



The proof of Theorem [T] is very similar to that of Theorem 1 in Glynn and Whitt (1992). 
A minor modification is necessary due to the relative nature of the stopping rule T^(e). 

Define z = zs/2 and V(n) = 2za n /y/n + p(n), where p{n) = o(ra -1 / 2 ). Then, 



13(e) = inf in > : 2za n / \fn + p(n) < eA n j 



can be denoted as 13(e) = inf jra > : V(n) < eA n j. Recall cr 2 . E (0, 00), then it is easy to 
verify that 

n 1 / 2 l/(n) -> 2zo e > w.p.l as n ->• 00. (8) 

By definition of 13(e), V{Tz(e) — 1) > eA T - 3 ( e )_i and there exits a random variable Z(e) E [0, 1] 
such that V(Ts(e) + Z(e)) < e\x 3 { t )+z(e)- Further note that T^(e) — > 00 w.p.l as e — > and 
hence A T3 ( e ) — > Ag w.p.l as e — > 0. Then using ^ we have 

limsupeT 3 (e) 1 / 2 < lim supTg^/V^e) - l)/A T3(e) _ 1 = 2za 6 /X e w.p.l. 



£->-0 

By a similar argument 



e-S>0 



lim inf eT 3 (e) 1 / 2 > lirr linf T^e) l l 2 V{T^{e) + Z(e))/A T3(e)+z(e) = 2za e /X e w.p.l. 



e^0 
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Thus, we have 

limeT 3 (e) 1/2 = 2za d /\ e w.p.l. (9) 



Given a FCLT at ([!]) holds and a n —¥ gq w.p.l as n — > oo, we have 

vV<7n (*n " On) 4 N(0, 1) (10) 
From ([9]) and ( 10 ) , it follows a standard random-time-change argument (p. 151 of Billingsley 



1995[ ) that 

V 7 ^)/^) (0 Ta (e) ~ 0^ 4 N(0, 1) as e -> 0. 

Finally, we have 



Pr (0 n e C[T 3 (e))) = Pr (9 Ta{e) - W G (-^^/^(e), ^ T3(f) / V^OO), 

= Pr (v / W/^T 3 (e)(^T 3 ( e ) - 0„)) € (-z, z)) 1 - <5 as e -)• 0. 
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